

%% "true" DGP point estimates
load('../true_DGP/tDGP.mat','DGP')
true_beta=DGP(1:(2*N));
clear DGP


%% reset rng seed
rng('default')


%%
NSP=1e+3; % number of simulated sample paths
T_dprk=4;
epsilon=mvnrnd(zeros(N,1),Sigma,(T-(54-T_dprk-1))*NSP)'; % growth shocks
epsilon=reshape(epsilon,N,T-(54-T_dprk-1),NSP);


%% small shock (6.5%) to Mongolia, Taiwan, India, Japan, South Korea, Philippines, (2004-T_chn):2004
P_dprk=zeros(size(P(:,:,54-T_dprk:T-1)));
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_dprk=y;
    y_dprk([131:132,134:136,148],54-T_dprk:54)=0.065+ep([131:132,134:136,148],1:T_dprk+1);
    PP_dprk=zeros(size(P(:,:,54-T_dprk:T-1))); %#ok<*PFBNS>
    BB_dprk=zeros(2*N,T-54+T_dprk);
    % 2001
    nm=find(M(:,54-T_dprk)==1);
    DD=[diag(1-D(nm,54-T_dprk)),diag(D(nm,54-T_dprk))];
    PP_dprk(:,:,1)=P(:,:,54-T_dprk);
    PP_dprk([nm;N+nm],[nm;N+nm],1)=PP_dprk([nm;N+nm],[nm;N+nm],1)+DD'*(Sigma(nm,nm)\DD);
    Pd_dprk=sqrt(diag(PP_dprk(:,:,1)\eye(2*N)));
    BB_dprk(:,1)=B(:,54-T_dprk);
    BB_dprk([nm;N+nm],1)=BB_dprk([nm;N+nm],1)+PP_dprk([nm;N+nm],[nm;N+nm],1)\(DD'*(Sigma(nm,nm)\(y_dprk(nm,54-T_dprk)-DD*BB_dprk([nm;N+nm],1))));
    U0_dprk=exp(alpha(133,1)+theta(1)*(Pd_dprk(133,1)*sgi_nodes(:,1)'+BB_dprk(133,1)+...
        ep_sd(133,1)*sgi_nodes(:,2)'));
    U0_dprk=(U0_dprk./(1+U0_dprk))*sgi_weights;
    U1_dprk=exp(alpha(133,1)+theta(2)*(Pd_dprk(N+133,1)*sgi_nodes(:,1)'+BB_dprk(N+133,1)+...
            ep_sd(133,1)*sgi_nodes(:,2)')-f(133,1));
    U1_dprk=(U1_dprk./(1+U1_dprk))*sgi_weights;
    DD_dprk=D;
    DD_dprk(133,54-T_dprk+1)=M(133,54-T_dprk+1).*(U1_dprk>=U0_dprk);
    for t=1:T-54+T_dprk-1
        if DD_dprk(133,54-T_dprk+t)~=D(133,54-T_dprk+t)
            y_dprk(133,54-T_dprk+t)=true_beta(133)*(1-DD_dprk(133,54-T_dprk+t))+true_beta(N+133)*DD_dprk(133,54-T_dprk+t)+ep(133,t+1);
        end
        BB_dprk(:,t+1)=BB_dprk(:,t);
        PP_dprk(:,:,t+1)=PP_dprk(:,:,t);
        nm=find(M(:,54-T_dprk+t)==1);
        DD=[diag(1-DD_dprk(nm,54-T_dprk+t)),diag(DD_dprk(nm,54-T_dprk+t))];
        PP_dprk([nm;N+nm],[nm;N+nm],t+1)=PP_dprk([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd_dprk=sqrt(diag(PP_dprk(:,:,t+1)\eye(2*N)));
        BB_dprk([nm;N+nm],t+1)=BB_dprk([nm;N+nm],t+1)+PP_dprk([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y_dprk(nm,54-T_dprk+t)-DD*BB_dprk([nm;N+nm],t))));
        U0_dprk=exp(alpha(133,1)+theta(1)*(Pd_dprk(133,1)*sgi_nodes(:,1)'+BB_dprk(133,t+1)+...
            ep_sd(133,1)*sgi_nodes(:,2)'));
        U0_dprk=(U0_dprk./(1+U0_dprk))*sgi_weights;
        U1_dprk=exp(alpha(133,1)+theta(2)*(Pd_dprk(N+133,1)*sgi_nodes(:,1)'+BB_dprk(N+133,t+1)+...
                ep_sd(133,1)*sgi_nodes(:,2)')-f(133,1));
        U1_dprk=(U1_dprk./(1+U1_dprk))*sgi_weights;
        DD_dprk(133,54-T_dprk+t+1)=M(133,54-T_dprk+t+1).*(U1_dprk>=U0_dprk);
    end
    D_dprk(:,:,nsp)=DD_dprk;
    B_dprk(:,:,nsp)=BB_dprk;
    P_dprk=P_dprk+PP_dprk/NSP;
end
D_dprk_ss=1*(mean(D_dprk,3)>0.5);
B_dprk_ss=mean(B_dprk,3);

% for Figure 9
Pinv_dprk=P_dprk;
for t=1:size(P_dprk,3)
    Pinv_dprk(:,:,t)=P_dprk(:,:,t)\eye(2*N);
end
unc_dprk=zeros(N,size(P_dprk,3));
for n=1:N
    for t=1:size(P_dprk,3)
        unc_dprk(n,t)=sqrt(Pinv_dprk(N+n,N+n,t)+Pinv_dprk(n,n,t)-2*Pinv_dprk(N+n,n,t));
    end
end
writetable(table(unc_dprk(133,:)),'../../Figures/Figure_9/dprk_unc_cfactual065.csv')


%% large shock (11%) to Mongolia, Taiwan, India, Japan, South Korea, Philippines, (2004-T_chn):2004
P_dprk=zeros(size(P(:,:,54-T_dprk:T-1)));
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_dprk=y;
    y_dprk([131:132,134:136,148],54-T_dprk:54)=0.11+ep([131:132,134:136,148],1:T_dprk+1);
    PP_dprk=zeros(size(P(:,:,54-T_dprk:T-1))); %#ok<*PFBNS>
    BB_dprk=zeros(2*N,T-54+T_dprk);
    % 2001
    nm=find(M(:,54-T_dprk)==1);
    DD=[diag(1-D(nm,54-T_dprk)),diag(D(nm,54-T_dprk))];
    PP_dprk(:,:,1)=P(:,:,54-T_dprk);
    PP_dprk([nm;N+nm],[nm;N+nm],1)=PP_dprk([nm;N+nm],[nm;N+nm],1)+DD'*(Sigma(nm,nm)\DD);
    Pd_dprk=sqrt(diag(PP_dprk(:,:,1)\eye(2*N)));
    BB_dprk(:,1)=B(:,54-T_dprk);
    BB_dprk([nm;N+nm],1)=BB_dprk([nm;N+nm],1)+PP_dprk([nm;N+nm],[nm;N+nm],1)\(DD'*(Sigma(nm,nm)\(y_dprk(nm,54-T_dprk)-DD*BB_dprk([nm;N+nm],1))));
    U0_dprk=exp(alpha(133,1)+theta(1)*(Pd_dprk(133,1)*sgi_nodes(:,1)'+BB_dprk(133,1)+...
        ep_sd(133,1)*sgi_nodes(:,2)'));
    U0_dprk=(U0_dprk./(1+U0_dprk))*sgi_weights;
    U1_dprk=exp(alpha(133,1)+theta(2)*(Pd_dprk(N+133,1)*sgi_nodes(:,1)'+BB_dprk(N+133,1)+...
            ep_sd(133,1)*sgi_nodes(:,2)')-f(133,1));
    U1_dprk=(U1_dprk./(1+U1_dprk))*sgi_weights;
    DD_dprk=D;
    DD_dprk(133,54-T_dprk+1)=M(133,54-T_dprk+1).*(U1_dprk>=U0_dprk);
    for t=1:T-54+T_dprk-1
        if DD_dprk(133,54-T_dprk+t)~=D(133,54-T_dprk+t)
            y_dprk(133,54-T_dprk+t)=true_beta(133)*(1-DD_dprk(133,54-T_dprk+t))+true_beta(N+133)*DD_dprk(133,54-T_dprk+t)+ep(133,t+1);
        end
        BB_dprk(:,t+1)=BB_dprk(:,t);
        PP_dprk(:,:,t+1)=PP_dprk(:,:,t);
        nm=find(M(:,54-T_dprk+t)==1);
        DD=[diag(1-DD_dprk(nm,54-T_dprk+t)),diag(DD_dprk(nm,54-T_dprk+t))];
        PP_dprk([nm;N+nm],[nm;N+nm],t+1)=PP_dprk([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd_dprk=sqrt(diag(PP_dprk(:,:,t+1)\eye(2*N)));
        BB_dprk([nm;N+nm],t+1)=BB_dprk([nm;N+nm],t+1)+PP_dprk([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y_dprk(nm,54-T_dprk+t)-DD*BB_dprk([nm;N+nm],t))));
        U0_dprk=exp(alpha(133,1)+theta(1)*(Pd_dprk(133,1)*sgi_nodes(:,1)'+BB_dprk(133,t+1)+...
            ep_sd(133,1)*sgi_nodes(:,2)'));
        U0_dprk=(U0_dprk./(1+U0_dprk))*sgi_weights;
        U1_dprk=exp(alpha(133,1)+theta(2)*(Pd_dprk(N+133,1)*sgi_nodes(:,1)'+BB_dprk(N+133,t+1)+...
                ep_sd(133,1)*sgi_nodes(:,2)')-f(133,1));
        U1_dprk=(U1_dprk./(1+U1_dprk))*sgi_weights;
        DD_dprk(133,54-T_dprk+t+1)=M(133,54-T_dprk+t+1).*(U1_dprk>=U0_dprk);
    end
    D_dprk(:,:,nsp)=DD_dprk;
    B_dprk(:,:,nsp)=BB_dprk;
    P_dprk=P_dprk+PP_dprk/NSP;
end
D_dprk_ls=1*(mean(D_dprk,3)>0.5);
B_dprk_ls=mean(B_dprk,3);


%% for Figure 9

% mean beliefs
writetable(table(B(N+133,:)-B(133,:)),'../../Figures/Figure_9/dprk_original.csv')
writetable(table(B_dprk_ss(N+133,:)-B_dprk_ss(133,:)),'../../Figures/Figure_9/dprk_cfactual065.csv')
writetable(table(B_dprk_ls(N+133,:)-B_dprk_ls(133,:)),'../../Figures/Figure_9/dprk_cfactual11.csv')

% belief uncertainty
writetable(table(uncertainty(133,:)),'../../Figures/Figure_9/dprk_unc_original.csv')
Pinv_dprk=P_dprk;
for t=1:size(P_dprk,3)
    Pinv_dprk(:,:,t)=P_dprk(:,:,t)\eye(2*N);
end
unc_dprk=zeros(N,size(P_dprk,3));
for n=1:N
    for t=1:size(P_dprk,3)
        unc_dprk(n,t)=sqrt(Pinv_dprk(N+n,N+n,t)+Pinv_dprk(n,n,t)-2*Pinv_dprk(N+n,n,t));
    end
end
writetable(table(unc_dprk(133,:)),'../../Figures/Figure_9/dprk_unc_cfactual11.csv')



